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ABSTRACT 

We explore the heating of the velocity distribution in the solar neighbourhood by 
stochastic spiral waves. Our investigation is based on direct numerical integration of 
initially circular test-particle orbits in the sheared sheet. We confirm the conclusion of 
other investigators that heating by spiral structure can explain the principal features of 
the age- velocity dispersion relation and other parameters of the velocity distribution in 
the solar neighbourhood. In addition, we find that heating by strong transient spirals 
can naturally explain the presence of small-scale structure in the velocity distribution 
("moving groups"). Heating by spiral structure also explains why the stars in a single 
velocity-space moving group have a wide range of ages, a result which is difficult to 
understand in the traditional model that these structures result from inhomogeneous 
star formation. Thus we suggest that old moving groups arise from irregularities in 
the Galactic potential rather than irregularities in the star-formation rate. 

Key words: solar neighbourhood - Galaxy: kinematics and dynamics - Galaxy: 
fundamental parameters - stars: kinematics - galaxies: kinematics and dynamics 
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The velocity distribution function (df) of stars in the solar 
neighbourhood provides unique insights into the Galactic 
potential field, the dynamical history of the Galactic disk, 
and the relationships between kinematics, age, and metal- 
licity for disk stars. 

Let us define the Local Standard of Rest (LSR) to be 
a fictitious point that coincides with the Sun at the present 
instant and travels in a circular orbit around the Galactic 
centre. We introduce a rotating Cartesian coordinate system 
with origin at the LSR, z-axis pointing radially outward, 
j/-axis pointing in the direction of Galactic rotation, and z- 
axis pointing to the south Galactic pole. The DF /(v) is de- 
fined so that /(v)dv is the number of stars per unit volume 
with velocities in the range [v, v + dv]. The st andard empir- 
ical m odel for the past century has been the ISchwarzschilcO 
il907l) DF, which is a triaxial Gaussian of the form 



/(v) oc exp 



(1) 



where v = (v x ,v y ,v z ). The two lowest moments of the DF 
are the mean velocity, v = (v x ,v y ,v z ), and the velocity- 
dispersion tensor, 
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where X = J f(v)X(v)dv/ J /(v)dv, and the matrix of^ is 
the inverse of the matrix oiij. 

In a steady-state, axisymmetric galaxy (i) the mean 
velocity relative to the LSR is tangential, so v x — v z — 
0; (ii) the tensors ay and of, are both diagonal; (iii) 
the ratio o~ xx /ay y is determined by the local gravitational 
force and its radial gra dient (e.g. IChandrasekharl 1 19421 . 
iBinnev fc Tremaindll987l see also eg. 1 141 below k 

The velocity dispersion of stars in the solar neigh- 
bourhood increases with age, probably because the disk is 
"heated" by one or more dynamical mechanisms. However, 
the interpretation of the observed age- velocity dispersion re- 
lation (AVR) is uncertain: (i) One school models the AVR as 
a smooth power law, a xx (t) oc t p , and interprets this be- 
havior as evidence for a continuous heating mechanism. Es- 
timates of the exponent p in the literature span the range 
0.2-0.5. (ii) Some authors argue that the dispersion rises 
steeply with age for stars < 5 Gyr old, and thereafter is rel- 
atively flat flCarlberg et al.lll985t iGomez et al.lll99^) . per- 
haps because the continuous heating mechanism saturates 
once the dispersion is large enough, (iii) A third model is 
that the dis persion does no t inc rease smoothly with age . 
For example. iFreemarJ dl99 if) and lQuillen fc Garnettl j200ll) 
argue that there are three discrete age groups (t < 3 Gyr, 
3 Gyr < t < 10 Gyr, t>10Gyr) with different dispersions, 
but within each group there is no evidence for a correla- 
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tion between dispersion and age. Such groups might arise 
because the continuous heating saturates after only 3 Gyr, 
and the higher dispersion of the oldest stars is due to a dis- 
crete event such as a merger. 

A wide variety of mechani sms for disk heating has 
been discussed (see lLacevI Il99ll for a review) : (i) Spitzer 
& Schwarzschild (1951, 1953) suggested that massive gas 
clouds could gravitationally scatter stars, leading to a steady 
increase in velocity dispersion with age, and thereby pre- 
dicted the existence of giant molecular clouds (GMCs). How- 
ever, heati ng by GMCs alone fails to explain several ob- 
serva ti ons llLace vlll984t IVillum senlll985| : I.Tenkins fc Binnevl 
Il990t lLacevI Il99lt IJenkinj Il992tl : the predicted ratio 
cr^/ffn of vertical to radial dispersion may be too high, 
roughly 0.72 compa red to the observed va l ue of 0.5 for 
old stars (but see llda. Kokubo fc Makinol Il993l for an 
opposing view); the predicted exponent in the AVR is 
somewhat too low, p<0.25; and the masses and num- 
ber density of GMCs determined from CO observations 
are too low to explain the observed dispersion, proba- 
bly by a factor of five or so but with substantial uncer- 
tainty, (ii) Transient spiral waves lead to potential fluctu- 
ations hi_th«_ohsk_tjiat_^c^ 

stars fearbanis fc Woltie3ll967t ISellwood fc Carlberdll984l 
ICarlberg fc Sellwoodlll985ir However, spiral waves only ex- 
cite the horizontal (x and y) velocity components effec- 
tively, since their characteristic spatial and temporal scales 
are much larger than the amplitude or period of verti- 
cal oscillations for young stellar populations, (iii) These 
considerations lead naturally to a hybrid model, in which 
spiral waves excite non-circular velocities in the plane, 
and the velocities are then redistributed between horizon- 
tal and vertical motion through GMC scattering iCarlberd 
Il987|). The hybrid model ha s bee n investigated thoroughly 
bv I.Tenkins fc Binnevl l|l99ff) and I.TenkinsI (jl992T ) using the 
Fokker-Planck equation. They find that they can success- 
fully reproduce the observed axis ratio a zz /a xx , the expo- 
nent p in the AVR and the radial dispersion of old stars, (iv) 
Other possible heating mechanisms, all of which rely to some 
extent on hypothetical or poorly understood components or 
features of the Galaxy, include scatterin g by massive com- 
pact halo objects or halo substruc ture llLacev fc Ostrikerl 
1985) , mergers with dwarf galaxi es jToth fc Ostrikerlll992t 
Iwalker. Mihos fc Hernauistll994lH~uang fc Carlberdll997t) . 
or the outer Lindblad resonance from the Galactic ba r 
jKalnaisHl99ll:lDehnenll999ll2000tlFuxl200ll:lQuilleJ2003i) . 

The Schwarzschild DF Q only approximates the veloc- 
ity distribution on the largest scales in velocity space. On 
smaller scales, there is substructure, which is most promi- 
nent in the youngest stars but present in stars of all ages. 
Discussion of substructure in the velocity DF dates back 
to Kapteyn's (1905) model of "two star streams", and for 
decades Eggen advocated the case for substructure in the 
form of "moving groups" in the solar neighbourhood llEggenl 
Il99fil and references therein) . 

Eggen and others have usually explained moving groups 
as the result of inhomogeneous star formation in the disk: 
in this model, stars in a moving group are born at a com- 
mon place and time, and then disperse into a stream that 
happens to intersect the solar neighbourhood. This model 
predicts that stars in a moving group should have the same 
age, metallicity, and azimuthal velocity (i.e. the same an- 



gular momentum, since this determines their mean angular 
velocity). The existence and membership of these groups 
was controversial until the Hipparcos satellite measured re- 
liable distances and proper motions for a large, homoge- 
neous database of nearby stars, and verified the presence of 
rich substructure in the velocity DF of both young and old 
stars, including a number of features that coin cide with mov- 
ing groups already identified by Eggen (e.g. |Pehner]|l998l 
Chereul et al. 1998,1999) 

In this paper we explore a quite different explanation 
for substructure in the velocity DF. We suggest that sub- 
structure arises naturally from the same spiral gravitational 
fluctuations that excite the growth of the velocity dispersion. 
In other words, substructure is caused by homogeneous star 
formation in an irregular potential, as opposed to inhomoge- 
neous star formation in a regular potential in the traditional 
model. 

We investigate this hypothesis by simulating the evolu- 
tion of the velocity DF induced by transient spiral structure 
in a simple two-dimensional model of the local Galaxy. We 
restrict ourselves to two dimensions because spiral structure 
does not excite vertical motions efficiently, and because the 
velocity DF ap pears to be well-mixed in the vertical direction 
l|DehnerJll998h . 

Following Eggen, we shall use the term "moving group" 
to denote substructure in the velocity DF of old stars ( > 1 
Gyr) at a given position. Unfortunately, the same term is 
sometimes also applied to OB associations, which are spa- 
tially localized conce ntrations of much younger stars (e.g., 
Ide Zeeuw et~ai]|l999Fl . 

Section[3|describes our simplified dynamical model. The 
results of our simulations are analyzed and compared to 
observed data in Section |S] examines briefly the tradi- 
tional hypothesis that substructure arises from inhomoge- 
neous star formation. Section [5] contains a brief discussion 
of the closely related process of radial migration of stars, 
and Sj7| contains concluding remarks. 



2 PROPERTIES OF SPIRAL STRUCTURE 

Our model depends on several properties of the Galaxy's spi- 
ral structure, such as the number of arms, the arm strength, 
and the pitch angle. We are concerned with spiral structure 
in the disk surface density (rather than, say, in the distribu- 
tion of young stars or gas). The properties of this structure 
in external galaxies are revealed by near-infrared images, 
which are do minated by the stars that c ontribute mos t of 
the mass jRix fc Riekdil993l but see also lRhoadsjll998l and 
Ijames fc _JMgar|ll999| f° r qu alifications) . 

iRix fc Zaritskvl l|l995h examined the K-band spiral 
structure of 18 face-on galaxies. They found that almost half 
had strong two-arm (m = 2) spirals, with arm-interarm con- 
trasts I m ax/-fmin — 2. If the arms are sinusoidal, with frac- 
tional amplitude e relative to the axisymmetric background, 
the arm-interarm contrast is l ma x/^min = (1 + e)/(l — e), so 
a co ntrast of 2 corresponds to e ~ 0.3. 

iBlock fc Pueraril Jl999T) examined K-band spiral struc- 
ture in 19 spirals. They found pitch angles ranging from 8° 
to 49° with a median value of 22°, and m = 2 amplitudes 
ranging from e = 0.03 to e = 0.5, with a median e = 0.1. 
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ISeiear fc James! (I1998T) conducted a similar survey of 
45 galaxies. They found that the dominant Fourier mode 
usually has m = 1 (36%) or m — 2 (31%), and the median 
pitch angle of the spiral structure was 8°, with little or no 
correlation with Hubble type. They measure the strength 
of the spiral arms in terms of the "equivalent angle"; their 
median equivalent angle of 14° corresponds to an amplitude 
e ~ 0.07 for a sinusoi d al m = 2 spiral. 

lElmeereen et all (Il999h stress that the amplitude of 
m — 2 spiral structure in the near-infrared depends on 
whether the o ptical spiral arms are classified as flocculent or 
grand-design felmeg rccn 1998). Grand-design spirals have 
arm-interarm contrasts of 1.5-6, corresponding to amplitude 
e = 0.2-0.7, while flocculent galaxies have contrast < 1-7, 
corresponding to e < 0.25. 

At optical wavebands, iMa et all l|l999l) find that the 
mean pitch angle for 51 Sbc galaxies (the same Hubble type 
as the Galaxy) is 15°. 

The measurement of spiral structure in our o wn Galaxy 
is mor e difficult than in face-on external galaxies. iDrimmell 
( 2000) uses K-band photometry of the Galactic plane to 
conclude that the Ga l axy c ontains a two-arm spiral with 
pitch angle 18° . IValiiel l|2002f l reviews a number of studies of 
the Galaxy's spiral structure, mostly based on young stars, 
gas and dust; these yield a range of pitch angles from 6°-20° , 
but Vallee concludes that the best overall fit is provided by 
an m — 4 spiral with pitch angle of 12°. The arm structure 
in the Galaxy appear s to be intermed iate between grand- 
design and flocculent l)Elmegreenlll998h . 



3 A NUMERICAL MODEL OF DISK HEATING 

We model disk heating in the sheared sheet, which approx- 
imates the local dynamics of a differentially rotating disk 
iSpitzer fc Schwarzschildl Il953l : iGoldreich fc Lvnden-Belil 
ll965l:IJulian fc Toomrell966l) . The LSR is assumed to travel 
around the Galactic centre in a circular orbit of radius Ro 
and angular frequency Q > 0. We use the same Cartesian 
coordinate system (x,y,z) introduced in the last Section, 
restrict ourselves to the z = plane, and denote position 
and velocity by x = (x,y) and v = (x,y). For \x\, \y\ <C Ro 
the equations of motion of a test particle are 



x - 2Q,y - ASlAx 
y + 20.x 



9$ 

dx ' 

dy ' 



(3) 



where $(x, t) is the gravitational potential due to sources 
other than the axisymmetric Galactic disk, and A > is the 
Oort constant 



A = 



2 \ n dRj Ro 



(4) 



where Rtt(R) is the circular speed at radius R and fi(-Ro) = 

n. 

If $ = the trajectories governed by equations have 
two integrals of motion related to energy and angular mo- 
mentum, 



E = \{x 2 + y 2 -AQAx 2 



H = y + 2Qx. 



(5) 



The solutions of the equations of motion are 
x = x g + a cos(k£ + (f>), 
V = Vg(t) asm(Kt + <j)), 

Vg(t) = y g o - 2Ax g t, 



(6) 



where [x g ,y g (t)] are the coordinates of the guiding centre, 
a is the epicycle amplitude, is a phase constant, and k is 
the radial or epicycle frequency, 



k = n 4 



or k = 4f2(f2 — A). 



(7) 



d In R J n 

The energy and angular momentum are related to the 
guiding-centre radius and epicycle amplitude by 



E = 2A(A-tyx 2 g + ±K 2 a 2 , H = 2(0.- A)x g 



(8) 



and the guiding-centre radius is related to the phase-space 
coordinates by 

y + 2Qx 



9 ~ 2(n-A) - 

The epicycle energy is defined as 



0) 



E x 



i r-2 , 
a l x + 



K ( 3) X n 



E + 



AH 



2{Q - A) 



12 2 
2 K a 

\± 2 + ^-(y + 2Ax) 2 . 



(10) 



The epicycle energy is closely related to the radial action 
/ = i K a 2 = E a /K. (11) 
For a particle in a circular orbit, 



(x,y,x,y) = {x g ,y g0 - 2Ax g t, 0, -2Ax g 



(12) 



thus E x = 0. 

The approximations used in deriving the linearized 
equations of motion © are only marginally valid in the solar 
neighbourhood, since the epicycle amplitudes a can be a sig- 
nificant fraction of Ro (for a population with radial velocity 
dispersion a xx in a galaxy with a flat rotation curve, 



&xx/Q 



(13) 



Thus old disk stars in the solar neighbourhood, with a xx — 
40kms _1 , have (a?) 1/2 ~ 0.18-Ro). Nevertheless, we believe 
that the sheared sheet accurately captures the most impor- 
tant features of the evolution of disk-star kinematics. 

The kinematics of a population of stars is described by 
the DF /(x, v, t), where /dxdv is the number of stars in the 
interv al (x, v) — > (x + rfx, v + rfv ). According to Jeans's the- 
orem feinnev fc Tremaindll987h . a stationary DF can only 
depend on the integrals of motion E and H (or E x , x g , a, 
etc.). In the solar neighbourhood [x = y = 0), the integrals 
are H = y and E x = \x 2 + 2Q 2 y 2 / k 2 . Thus in a steady state 
the velocity distribution must be an even function of x. 

It is straightforward to show that the mean velocity and 
velocity-dispersion tensor of any stationary, spatially homo- 
geneous DF in the sheared sheet must satisfy the relations 



0, Vy 



*2Ax, <7 X y 



0. 



n-A 



(14) 
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A useful model DF for the sheared sheet is 

/(x,v) oc exp (-E x /<7o) , (15) 

where E x is defined in equation HUH . At a given location, 
the DF 115t leads to a two-dimensional version of the triaxial 
Gaussian velocity distribution JQ, in which the mean veloc- 
ity and velocity-dispersion tensor satisfy the relations (1141) . 
Thus the DF 11511 is also sometimes (confusingly) called the 
Schwarzschild DF. 

In our simulations the relations 114^ are not satisfied 
exactly, so it is useful to introduce the principal-axis system 
(xi,X2) in which the velocity-dispersion tensor is diagonal; 
the xi-axis is chosen to lie within 45° of the i-axis. The 
corresponding velocity dispersions are o\ and 02, which we 
will normally plot instead of a xx and a yy \ usually in our 
simulations and in the data we shall find that ai > 02. The 
vertex deviation l v is the Galactic longitude of the xi-axis, 
and is given by 



lv 



■ arctan , 2 

« \ 0° XX 



2ol 



(16) 



where |Z„| < 45°. Stellar populations in the solar neigh- 
bourhood have vertex dev iations in the range < l v < 30° 
jBinnev fc Merrifieldll99gl) . 

Throughout the paper, we assume that the rota- 
tional curve of the underlying axisymmetric galaxy is flat, 
Rfl(R) = constant, so that A = if2 and n = \J~2£l. We 
shall assume that the surface density of the disk in the 
solar neighbourhood is Ed = 50Mq pc ~ 2 ; recent observa - 
tional estimates are S d = 40M p pc~ 2 JCreze et allll998fr . 
48M w pc~ 2 faolmberg fc Flvnnll200ol) . and 42 ± 6M pc -2 
iKorchagin et alTl2003T) ~ The assumed surface density only 
enters our analysis through the definition of the fractional 
spiral amplitude e below (eq. l21ll . 



3.1 Spiral waves 

We approximate spiral structure as a superposition of waves 
with surface density and potential 



T,(x,y,t) = E s (i) exp i(k x x + k y y), 
$(x,y,t) = & s (t) expi(k x x + k y y), 



(17) 



where k = (k x ,k y ) is the wavenumber, and only the real 
part of E or $ is physical. Without loss of generality we can 
assume that k y > 0; then the spiral is trailing if k x > and 
leading if k x < 0. In this paper, we only consider trailing 
spiral waves. The number of arms m and the pitch angle a 
are determined by k through the relations 



111 



tan a. 



(18) 



The relation between p otential and surface densit y for spiral 
waves is given by (e.g. iBinnev fc Tremaindll987h 



27rGS 3 (t) 
k ' 



(19) 



where k = |k| . 

Steady spirals, in which E s (t), <E> s (t) oc exp(iu}t), heat 
stars on nearly circular orbits only a t the Lindblad reso- 
nances jLvnden- Bell fc Kalna'i3ll972D . which occur at the 
guiding centre radii given by 



(20) 

We focus instead on transient spiral structure, which can 
heat stars over a range of radii. We model each transient 
spiral using a Gaussian amplitude dependence centred at 
time t s with standard deviation a 3 : 



*.(*) 



2neGT, d 
k 



■ exp 



2al 



+ i(6 + 2Ak y x c t) 



(21) 



Here 6 is a phase constant, and x c is the corotation radius 
relative to the LSR (the phase of the spiral is constant for an 
observer on a circular orbit at x c ). The parameter e measures 
the amplitude of the spiral, and the normalizing constants 
are chosen so that the maximum surface density in the spiral 
is a fraction e of the surface density of the underlying disk 
(eci.im 

In each simulation the trajectories of the stars were fol- 
lowed between time t — and t = to. During this interval 
the stars were perturbed by N„ transient spiral waves. Each 
wave had phase constant 6 chosen randomly from [0, 27r], 
corotation radius x c chosen randomly from a Gaussian dis- 
tribution centred on the LSR with standard deviation a c , 
and central time t s chosen randomly from [— 2a s ,to + 2a 3 ] 
(this is slightly longer than the integration interval, to in- 
clude the effects of transients whose wings are inside the 
integration interval although their central times are not). 

The fractional amplitude of the surface-density per- 
turbation due to a single transient spiral is e at the wave 
peak. However, in some ways a better quantity to compare 
with the observational data on spiral amplitudes is the root- 
mean-square (RMS) time average of the fractional surface- 
density amplitude, 



ri/2 



N a a s 



to 



1/2 



(22) 



A closely related quantity is the RMS potential perturbation, 



27rGe r msE ( ii?0' 



V2r, 



(23) 



(the factor %/2 arises because <Ev ms is the RMS fluctua- 
tion rather than the RMS am plitude, which is smaller by 
V2). Jenkins & Binncv (1990) use Fokker-Planck calcula- 
tions of heating by transient spiral structure to estimate 
that $ rms = (9-13kms -1 ) 2 . 

The power spectrum of each transient is a Gaussian 
with standard deviation (2o- 2 )~ 1//2 . The central frequencies 
of the power spectra of the transients follow a Gaussian dis- 
tribution with standard deviation 2Ak y a c . Thus the power 
spectrum of the ensemble of N a transients is smooth if the 
overlap factor 

-2\-l/2 



C = N s 



2ai 



2 3 / 2 Ak y a c a a 



(24) 



is large compared to unity; on the other hand if C <C 1 the 
power spectrum and hence the heating is localized at narrow 
resonances. All of our simulations have C> 1 (see Tablc0. 

The dispersion of the power spectrum of all the tran- 
sients is at, given by 

a * = 2^ + ( 2Ak y a ^ 2 - ( 25 ) 

Note that in this model problem different azimuthal 
wavenumbers m and m' = fm axe equivalent if we rescale 
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the other variables by: k' = /k, a' = a, t' = t, (a/, y') = 
(x/f,y/f), $' s = $ s // 2 ,^ = E.//. e ' = < = 

(T c = Co//, (To = (To//, (Ts = (T s . 



3.2 Determining the solar neighbourhood velocity 
distribution 

Assume that the disk stars are formed on circular orbits at 
time 0. We wish to determine the velocity distribution at the 
present time to in the solar neighbourhood, (x,y) = (0,0). 
This is a two-point boundary-value problem rather than an 
initial-value problem: the boundary conditions on the phase- 
space coordinates are (x,y,x,y) = (x g , y g o, 0, — 2Ax g ) at 
t = (on circular orbits) and (x, y) = (0, 0) at t = to (in 
the solar neighbourhood). The solutions to the boundary- 
value problem will be a set of distinct initial positions 
[xi(t — 0),j/j(t = 0)] or final velocities [xi(t = to),Vi(t = 
to)] = (— Ui, Vi), i = 1, 2, . . . (the sign of it is chosen to agree 
with the usual convention that stars moving towards the 
Galactic centre have u > 0). Because the motion of stars in 
the fluctuating gravitational field of the transient spirals is 
complicated, we expect that there will be many solutions. 
In this idealized model the velocity-space DF in the solar 
neighbourhood will consist of a set of delta-functions at the 
velocities (m,Vi), but in practice these will be smeared into 
a continuous distribution by observational errors, the non- 
zero volume of the "solar neighbourhood" , the small initial 
velocity dispersion of the stars when they are formed, etc. 

Normally, boundary- value problems are solved most ef- 
ficiently by iterative methods. In this case, however, most 
of the solutions have only a tiny domain of attraction in 
the (u, v) plane; thus iterative methods are inefficient. We 
have therefore a dopted a differ ent approach, based on Monte 
Carlo sampling iDehnenll2000F) . 

Stars are not born on precisely circular orbits, in part 
because star-forming clouds are not on circular orbits. 
We may therefore assume that the stars initially have a 
Schwarzschild DF (eq. 1151 , with a small but non-zero ini- 
tial velocity dispersion o~o- 

Following Dehnen, we integrate orbits backward in 
time, starting at t = to (which we call the initial integration 
time or "present epoch") and ending at t = to (the final 
integration time or "formation epoch"). The initial condi- 
tions are chosen at random from (x,y,x,y) = (0,0,— u,v), 
with |m|, \v\ < Umax- We choose iWx = llOkms -1 , large 
enough to include almost all disk stars in the solar neigh- 
bourhood. We then integrate equations (J^J backwards in 
time to t he formation epoch. The collisionless Boltzmann 
equation (iBinnev fc Tremaindll987l) states that the phase- 
space density around a trajectory is time-invariant. There- 
fore the phase-space density Fi at (ui,Vi) at the present 
epoch is given by 1151 . where the epicycle energy E x is 
measured at the formation epoch. We assume that the spa- 
tial volume of our survey of the solar neighbourhood is in- 
dependent of velocity, so the density in velocity space at 
the present epoch is proportional to Fi. Dehnen's procedure 
therefore provides a Monte Carlo sampling of the velocity- 
space DF at the present epoch. 

We convert our Monte Carlo realization to a smooth 
DF by convolution with observational errors. We replace the 
point solutions = (ui,Vi) by Gaussians, that is, 



/(v) = > RS(v-Vi 



is replaced by 



/(v) = 



exp 



( v ~ Vi. 



(26) 



(27) 



where <r b is the assumed observational error. 

We have described the procedure for estimating the DF 
of stars of age to, using equation 115H and the epicycle en- 
ergies at t — 0. However, we can also determine the DF for 
stars having an intermediate age to — t m in the course of 
the same orbit integration, using the epicycle energies of the 
stars at the intermediate time t m . This approach measures 
the dependence of the DF on stellar age, and thus determines 
the AVR for a given realization of the set of transient spiral 
arms, using only a single set of orbit integrations. 

3.3 Units and model parameters 

The time unit is chosen to be fi -1 and the distance unit 
is chosen to be Ro. The azimuthal orbital period is then 
2n/Q = 2tt. For conversion to physical units we assume that 
Ro = 8kpc and flRo = 220 km s -1 , so that the unit of time 
is 35.6 Myr. Most of our numerical integrations ran for an 
interval to = 281, corresponding to 40.8 rotation periods of 
the LSR or 10.0 Gyr. 

We consider only trailing spirals with m = 2 or m = 4. 
The standard deviation of the distribution of corotation 
radii for the transient spirals was <r c = 0.25 correspond- 
ing to 2 kpc. With these choices the overlap factor is 
C = \AlN s {m/2)/{Q.a s ). 

The radial velocity dispersion of stars at the time of 
their birth is taken to be ao = Skms^ 1 . The observational 
error in the velocities is taken to be a b = 3kms _1 ; for the 
typical Hipparcos parallax error of 1 mas, a ^ corresponds 
to the velocity error arising from the distance uncertainty 
for a star with transverse velocity 30kms _1 at a distance of 
100 pc. 

In each simulation we launch either 1 x 10 r or 7 x 10 7 
stars from the solar neighbourhood and integrate back for 
10 Gyr (the smaller simulations are used to estimate mo- 
ments, such as those in Figure Q and the larger simulations 
are used to estimate the DF, as in Figure [7|. In the larger 
simulations, typically 10000-15000 stars contribute signifi- 
cantly to the solution, in the sense that their epicycle energy 
at the formation epoch is < 3(Tq- This number is larger than 
the number of stars in the Hipparcos samples to which we 
compar e (e.g. 46 00 in the largest color-selected sample anal- 
ysed bv lDehnenll998f) . We have verified that the conclusions 
below are unaffected if we double the number of stars in the 
simulation or change the random-number seed used to se- 
lect the sample of stars (keeping the same seed to select the 
spiral transients). 

We carried out a number of simulations with different 
values of the spiral- wave parameters m (number of arms), N B 
(number of spirals), e (fractional surface-density perturba- 
tion), a s (duration of the transient), a c (dispersion in coro- 
tation radius), and a (pitch angle). The parameters of the 
simulations are shown in TableQ Each simulation is defined 
by these six parameters and the random number seeds for 
both the stars and the spiral arms. 
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Table 1. Simulation parameters 
Parameter Symbol Run 5 Run 10 Run 20 Run 40 Run m4 Run sim 



Defining parameters 



Number of spiral arms 


m 


2 


2 


2 


2 


4 


2 


Number of spiral waves 


N 3 


45 


18 


12 


9 


12 


48 


Maximum fractional surface 
















density of spiral waves 


c 


1.61 


0.80 


0.80 


0.67 


0.80 


0.31 


Time scale of spiral waves 


<Ts 


1 


1 


1 


1 


1 


1 


RMS corotation radius (R. — Ro)/Ro 




0.25 


0.25 


0.25 


0.25 


0.25 


0.25 


Pitch angle 


a 


5° 


10° 


20° 


40° 


20° 


20° 




Derived parameters 










Overlap factor (eg. 1241 


C 


64 


25 


17 


13 


34 


68 


RMS surface-density amplitude (eq. 1221 




0.86 


0.27 


0.22 


0.16 


0.22 


0.17 


RMS potential in (kms -1 ) 2 (eg. 1251 


^rms 


(16.9) 2 


(13.4) 2 


(17.0) 2 


(19.8) 2 


(12.0) 2 


(15.0) 2 



Our runs have pitch angles of 5°, 10°, 20°, and 40°, 
which span the range of observed pitch angles in spiral galax- 
ies. We then choose the amplitude e and the number of tran- 
sients N s to approximately reproduce the observed radial 
velocity dispersion of old main-sequence stars in the solar 
neighbourhood . 

All of our runs have two-armed spirals (m = 2), except 
for run rh4, which has the same parameters as run 20 except 
that m — 4. 

We also have run sim, in which we have used a larger 
number of weaker transients than in run 20, and have also 
adjusted some of the central times of the spiral transients 
to better reproduce the features of the observed solar neigh- 
bourhood DF (see 

The RMS potential fluctuation <f> rms = (13-20 km s' 1 ) 2 
in our four runs, highe r than the velu e $rm s = 
(9-13kms -1 ) 2 derived bv I.Tenkins fc Binnevl jlQQOT l from 
Fokker-Planck calculations by about a factor of 2-2.4. 
A minor component o f this difference arises because 
I.Tenkins fc Binnevl ([[990) include GMCs, which contribute 
about 20% of their total heating rate. A more important 
effect, pointed out to us by A. J. Kalnajs, is that in our 
model the heating is spatially inhomogeneous: our assump- 
tions place the Sun fairly close to corotation for most of 
the transients, so that heating at the Lindblad resonances is 
more effective inside and outside the solar circle than it is 
at the solar circle. 

The fractional surface-density fluctuation e required to 
reproduce the AVR in run 5 exceeds unity, and the RMS 
surface-density fluctuation is close to unity; these high values 
are unrealistic and suggest that a larger number of weaker 
transients is required if the pich angle is as small as 5°. 

Runs 10, 20, 40, and m4 have peak amplitudes e = 0.7- 
0.8, similar to the strongest grand-design spirals, but these 
large amplitudes are only present for a small fraction of the 
time (N a a 3 /to = 0.03-0.06). The relatively small number of 
high- amplitude transients in these runs is consistent with a 



1 More precisely, main-sequence stars redward of the Parenago 
discontinuity at B — V = 0.6 are all believed to have the 
same mean a ge, and have a color-ind ependent radial dispersion 
~ SSkrns" 1 foehnen fc BinnevllT99sl) . We match this to the ra- 
dial dispersion in our simulation by assuming a constant star- 
formation rate, which is consist ent with solar neighbourhood dat a 
within the large uncertainties feinnev. Dehnen fc BertellftOOd) . 



model in which the Galaxy has experienced several short- 
lived grand-design spiral phases, possibly caused by minor 
mergers. In contrast, run sim has a much larger number of 
weaker transients (e = 0.3, N s a s /to = 0.17). 

Each run was repeated seven times (realizations a-g), 
using different random-number seeds to generate both the 
stars and the spiral transients, in order to distinguish sys- 
tematic from stochastic variations. To study extremely long- 
term behavior, we also ran a realization h, in which both the 
integration time and the number of transients are a factor 
of three larger than in a-g. 



4 RESULTS 

4.1 Age- velocity dispersion relation 

In H3.2I we have shown that a single set of orbit integrations 
can be used to derive the DF for stars of all ages between 
and to. Thus we may combine DFs of different ages to pro- 
duce the DF that corresponds to a uniform distribution of 
ages between and 10 Gyr (hereafter the "combined DF"), 
as would result from a uniform star-formation rate. The av- 
erage moments of this DF for our runs are s hown in Tabled 
along with the observed values taken from iDehnenl |l998). 
Of course, the close agreement of the radial velocity disper- 
sion a xx in each run with the observed value of 37kms~ 
arises because the strength and number of spiral transients 
were chosen to produce the correct radial dispersion. 

The evolution of the principal axes ai and 02 of the 
velocity-dispersion tensor in runs 5, 10, 20 and 40 is shown in 
Figure [T] The three rows of panels show the results for three 
different realizations, i.e., three different sets of random- 
number seeds used to generate both the stars and the spiral 
transients. In general, the results in this and other figures 
are independent of the realization of the stellar distribution, 
and all of the variations seen are due to changes in the real- 
ization of the transient spiral structure — a reflection of the 
fact that the Monte Carlo simulation of the stellar distri- 
bution is a numerical method while the simulation of the 
distribution of spiral transients is a model of a stochastic 
physical process. 

The AVRs illustrate the following: 

• Open spirals (large pitch angle) heat the disk stars more 
effectively than tightly wound spirals. For example, both the 
number N 3 and the surface-density amplitude e of the spiral 
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Figure 1. The age dependence of the principal axes of the velocity-dispersion tensor in runs 5, 10, 20 and 40. The solid line shows a\ 
and the dashed line shows <T2- The three rows (a,b,c) represent different realizations. The dot-dashed lines show the best fit of equation 
1291 for each realization; see also Table 131 



Table 2. Simulation results and observations 
Quantity Symbol Run 5 Run 10 Run 20 Run 40 Run m4 Run sim Observations 



Radial velocity dispersion 


iTii(kms x ) 


36 ±1 


38 ±2 


40 ±4 


39 ±6 


35 ±3 


40 ±3 


36.9 


vertex deviation 


lv(°) 


-3 ±4 


3±3 


8 ±16 


-8 ± 12 


3 ± 8 


-6 ± 11 


10.9 


axis ratio of velocity ellipsoid 


0"l/<72 


1.3 ±0.1 


1.4 ±0.1 


1.2 ±0.1 


1.3 ±0.1 


1.4 ±0.1 


1.4 ±0.2 


1.48 


mean radial velocity 


u( kms -1 ) 


0±4 


-1±3 


0±5 


3±1 


0±4 


0±8 


-0.5 


mean azimuthal velocity 


v{ kms -1 ) 


0±2 


0± 2 


-1±3 


3±1 


-2 ±3 


2 ± 5 


-17.7 



Results given are the mean and standard deviation from realizations a— g, for a population o f stars with ages uniformly distributed between 
and 10 Gyr. "Observations" column gives results from samples B4+GI from Dolmen ( 1998), which contains stars redward of the Parenago 
discontinuity (B — V> 0.6, for which the lifetime of main-sequence stars exceeds the age of the Galaxy) and giant stars. The large negative 
mean azimuthal velocity in the observations is due to the asymmetric drift, which is not present in our simulations. 
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transients are larger in 5 than the other runs shown (see 
Table , even though all runs produce the same velocity 
dispersion by design. 

• The heating occurs in discrete steps separated by peri- 
ods of zero growth. This feature arises because we have only 
a small number of spiral transients over the course of the 
simulation, each of which has the same amplitude and lasts 
only for a characteristic time a s = 1 out of the total inte- 
gration time of 281. During the intervals between the spiral 
transients, the DF is stationary so the velocity dispersions 
are constant. 

• In run 5 the growth of the dispersion slows sharply after 
the radial dispersion reaches ~ 30kms _1 . This effect prob- 
ably arises because the epicycle amplitude becomes compa- 
rable to the radial spacing between the spiral arms, so that 
the forces from the spirals tend to average to zero over an 
epicycle oscillation. More specifically, the radial spacing of 
the spiral arms is Aa; = 2irRo tan a/m; the RMS epicycle 
amplitude (a 2 ) 1//2 is given by equation (1131 . and we expect 
the spiral heating to become ineffective when the RMS epicy- 
cle amplitude exceeds the peak-to-trough distance between 
arms, that is, when (a 2 ) 1 ^ 2 /Ax > 0.5. For small pitch angle 
a, 



2\1/ 2 



ma xx 



Ax 



2ttQ.Ro tan a 



0.50- 



(Tl 



5° 



2 SOkms" 1 a 



(28) 



• There are substantial differences among the realizations 
of a given run; in other words there are significant stochastic 
effects in spiral-wave heating (see the discussion of the ex- 
ponent p below for further detail). The stochastic variations 
in the heating rate become smaller as the velocity dispersion 
increases (see also Figure which follows the AVRs over an 
even longer time interval). 

For comparisons with theoretical models, it is useful 
to fit the simulated AVRs wit h a paramet rized function. A 
common parametrization (e.g. lLacevlll991l) is 



<n(t) = K 1/p + Ct) p , 



(29) 



where oo is the initial dispersion along the major principal 
axis, fixed at 3kms _1 as in H3.3I The motivation for this 
form is the assumption that the evolution of the dispersion 
is governed by a diffusion equation with velocity-dependent 
diffusion coefficient, 



(30) 



da 2 (t) 
dt 

where C\ = 2pC and q = p" 1 — 2, with C, p and q constants. 
For o\ (To, equation (12911 implies o\(t) oc t p . 

The dashed lines in Figure show the results of least- 
squares fits for the parameters C and p in equation (1291 . 
In each case, the fitting function successfully preserves the 
overall shape of the AVR but smooths over fluctuations due 
to the finite number of transients. We fit each of the seven 
realizations of each run and obtained the values of the AVR 
exponent p given in Tabled] The standard deviation in p can 
be as large as 35% (in run 20) among realizations of a single 
parameter set, and the variation due to changes in the spiral- 
structure parameters is just as large. Thus, the exponent of 
the AVR is not likely to discriminate well between heating by 
transient spiral structure and heating by other mechanisms. 

Figure [5] shows the AVRs of realization h, which was 



run for to — 30 Gyr, three times as long as the other 
runs. The purpose of this long run was to see how well the 
parametrization 1291 worked as the heating process contin- 
ued. The dashed lines show that equation 1291 still provides 
a reasonable fit; however, the best-fit exponent p (see Table 
|3J| is systematically smaller than in the shorter realizations 
a-g, except in run 5. Therefore, equation (1291 should be con- 
sidered as a fitting function valid over a limited time interval, 
rather than a formula with predictive power. 

All of our runs have the same value of the parameter 
<j s (Table 0, which represents the duration of the spiral 
transients. We have experimented with a range of values of 
a 3 . As shown in Figure [3] the dispersion for the combined 
DF with a uniform distribution of ages between and 10 
Gyr is almost independent of a s so long as a s > 0.5. A prob- 
able explanation is that the width of the power spectrum 
of the transient perturbations, given by equation 1251 . is 
dominated by the dispersion in pattern speeds rather than 
the duration of a single transient once <r s > l/(2V2Ak y a c ), 
which occurs when cr s > 1.4 for our parameters. Of course, 
when a s becomes very large, the overlap factor of equation 
1241 . will decrease below unity and the heating will be lo- 
calized at discrete resonances. 



4.2 The velocity ellipsoid 

The time evolution of the vertex deviation is shown in Fig- 
ure 2] Not surprisingly, the stochastic variations in the ver- 
tex deviation are larger for runs 10, 20 and 40, which have 
larger pitch angle and fewer transients, and the fluctuations 
in the vertex deviation decline with time as the growth rate 
of the velocity dispersion declines. The small but non-zero 
vertex deviations seen near the end of the runs are roughly 
consistent with the vertex deviation ~ 10° seen in the old 
stellar population in the solar neighbourhood. There is no 
evidence of a systematic preference for positive or negative 
vertex deviations for old stars. 

Figure shows the evolution of the ratio of the ma- 
jor and minor principal axes of the velocity ellipsoid. As in 
the case of the vertex deviation, the fluctuations are much 
larger in the runs with larger pitch angle. In a steady-state 
axisymmetric galaxy with a flat rotation curve, the axis ra- 
tio should be ai/a 2 = [0/(0 - A)] 1/2 = 2 1/2 = 1.414 (eq. 
1141 . and the observed axis ratio gradually approaches this 
value as the runs progress. 

Figure |5] shows the evolution of the mean radial and 
azimuthal velocity in the solar neighbourhood. There are 
significant fluctuations, up to lOkms -1 for the larger pitch 
angles, even in the oldest stellar populations. These fluc- 
tuations limit the accuracy with which we can construct 
axisymmetric models of local Galactic kinematics. 

We have also plotted the axis ratio of the velocity ellip- 
soid against the vertex deviation but have found no system- 
atic correlation between these two quantities. 

The parameters of the DF for a population of stars with 
a uniform distribution of ages between and 10 Gyr are 
shown in Table |5] They are all consistent with the observa- 
tions, except for the mean azimuthal velocity v. The discrep- 
ancy in v reflects the fact that our simulations do not show 
asymmetric drift because of the approximations inherent in 
the sheared sheet [eqs. (J^J with the perturbing potential 
$ = are symmetric under reflection through the origin]. 
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Figure 3. The dependence of the velocity dispersion on the duration of the transients, derived from the combined DF for a population 
of stars with a uniform distribution of ages up to 10 Gyr. All of the runs plotted in a given panel use the same defining parameters as 
in Table 111 except for the transient duration a s . The solid and dashed lines show a xx and cr yy , respectively. The velocity dispersion is 
almost independent of the transient duration a s , except for cr s <0.5. 



4.3 Shape of the velocity distribution 

Figure [7| shows the DF for stars with age 5 Gyr, from runs 
5a, 10a, 20a and 40a. We present the DF in two forms: (i) 
we display all of the points from our sampling of velocity 
space at the present epoch that yield nearly circular orbits at 
the formation epoch (i.e. those with epicycle energy < 3<tq, 
(Jo = 3kms _1 ). (ii) We show contour plots of the smooth 
DF that is expected with Hipparcos observational errors at ~ 



100 pc (see eg . 12711 . for two different realizations of the stellar 
distribution and the same realization of the spiral transients. 
The difference between these two realizations (the second 
and third rows) is an indication of the numerical noise in 
our simulations, which is small. 

The simulated DFs in the second and third row of pan- 
els do not resemble the Gaussian or Schwarzschild DF (eq. 
0. They are much lumpier, and the lumps are distributed 
more-or-less uniformly over a large area in velocity space, 
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Table 3. Fitted exponent p of AVRs Ceo. 1551. 





a 


b 


c 


d 


c 


f 


g 


a-g a 


h 


5 


0.31 


0.20 


0.22 


0.27 


0.27 


0.22 


0.31 


0.25±0.05 


0.24 


10 


0.30 


0.34 


0.37 


0.32 


0.24 


0.24 


0.25 


0.29±0.05 


0.25 


20 


0.33 


0.46 


0.35 


0.59 


0.29 


0.76 


0.57 


0.48±0.17 


0.29 


10 


0.37 


0.50 


0.38 


0.53 


0.36 


0.51 


0.40 


0.44±0.08 


0.23 


in 1 


0.48 


0.28 


0.26 


0.44 


0.24 


0.53 


0.26 


0.36±0.12 




sim 


0.40 


0.43 


0.55 


0.40 


0.39 


0.41 


0.66 


0.46±0.10 





The mean and standard deviation of the values of p from runs a- 
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Figure 4. The evolution of the vertex deviation l v (ea. 1161 . The three rows represent different realizations. 



rather than being concentrated at zero velocity. Qualita- 
tively, the Schwarzschild DFs resemble a contour map of a 
single large mountain, whereas the simulated DFs resemble 
an entire mountain range. The "lumpiness" or degree of sub- 
structure in the DF is a strong function of the pitch angle 
of the spiral arms, and in the range of pitch angles 10°-20° 



that observations suggest for our Galaxy the DFs are highly 
irregular. 

The DFs in Figure [7] exhibit prominent streaks in the 
direction v = v y =constant, particularly at the larger pitch 
angles. The explan ation of t hese s treaks is straightforward, 
and dates back to IWoollevI (Il96ll) . It is straightforward to 
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Figure 5. The evolution of the ratio o\/<J2 between the major and minor principal axes of the velocity-dispersion tensor. In a steady-state 
axisymmetric galaxy with a flat rotation curve, the axis ratio should be 2 1 / 2 = 1.414. 



show from equations © that for stars in the solar neigh- 
bourhood [x = y = 0) the radial and tangential velocities 
(it, v) are related to the guiding-center coordinates (x g ,y g ) 
by 



(x g ,y g ) = 



1 



2(n - A) 



{v,u). 



(31) 



Suppose that scattering from a spiral transient leads to 
an enhancement in the density of stars over the range of 
guiding-centre radii x g — * x g + Ax g . Since the guiding-centre 
velocity is y g = — 2Ax g (eq.HjJ, the stars will slowly spread 
out into an extended streamer in physical space; after time 
t the streamer will have width of order the epicycle ampli- 
tude and length Ay ~ 2AAx g t. The stars found at a given 
location along the streamer (for example the solar neigh- 
bourhood) will thus have very similar values of x g and thus 
of v. In other words, a tidal streamer will present itself in 
the local velocity distribution as a streak with constant v. 



The 1:1 relation between (u,v) and (y g ,x g ) implied by 
equation 1311 implies that the structure seen in Figure |7] is 
clo sely related to the guidin g center correlations discussed 
bv lToomre fe Kalnaisl il99ll) . 



ISkulian. Hearnshaw fc Cottrelll jl999T) find that the ob- 
served DF in the solar neighbourhood exhibits parallel 
"branches" of enhanced density in velocity space. It is 
tempting to id entify these with the streaks predicted by 
IWoollevI (1196 ll) and seen in Fig ure |7| Howe ver, th e branches 
found bv ISkulian. Hearnshaw fe Cottrelll ||l999j) are tilted 
with respect to the expected orientation v ^constant; they 
find instead that v = constant — 0.47m. Thus the correct 
interpretation of these branches is unclear: they may arise 
from some physical effect we have not modeled, or they may 
be artifacts arising from our predisposition to find order 
in random patterns. Note, for example, that parallel tilted 
"branches" appear to be present in the contour plots of the 
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Figure 6. The evolution of the mean velocities u = —v x (solid line) and v = v y (dotted line). 



DF in run 20a shown in Figure [7| with a slope about half of 
that claimed by Skuljan et al. 



tion is that they focused on a single steady bar-like potential 
whereas we have studied stochastic spiral transients. 



The principal conclusion from Figurc|7]is that the lumps 
in the simulated DFs, at least in the simulations with pitch 
angles > 10°, bear a striking resemblance to the moving 
groups in the solar neighbourhood (see for example iDehnenl 
ll99alSkulian. Hearnshaw fc Cottrdll999l or Chereul et al. 
1998, 1999), which strongly suggests that at least the older 
moving groups are a consequence of the same transient spiral 
structure that heats the stellar distribution. Thus we sug- 
gest that old moving groups arise from fluctuations in the 
gravitational potential rather than fluctuations in the star- 
formation rate, as has usually b een assumed. T his result was 
anticipated to some extent bv iKalnaisI lll99ll) and Dehnen 
(1999, 2000), who argued that some of the most prominent 
substructure in the solar neighbourhood DF arose from the 
outer Lindblad resonance with the Galactic bar. The distinc- 



4.4 Evolution of moving groups 

One of the striking features of the solar neighbourhood DF 
is that the same moving groups appear to be pr esent in 
stella r populations of different ages (see Table 2 of lDehnenl 
1998 or the top row of Figure |5J. This result is difficult to 
understand in the traditional view that moving groups arise 
from inhomogeneous star formation. Thus it is worthwhile 
to explore the temporal evolution of the DF and the survival 
of moving groups in our simulations. 

Our standard four simulations, 5, 10, 20, and 40, are not 
ideal for this task. As we have seen, run 5 has relatively weak 
moving groups and unrealistic parameters. The other runs 
show a strong moving group at the origin in the combined DF 
which is not present in the observed DF. This moving group 
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Figure 7. The velocity DF for stars of age 5 Gyr, from runs 5a, 10a, 20a, and 40a. The top row shows locations from our Monte 
Carlo sampling of velocity space that yield nearly circular orbits (epicycle energy < 3oq, r?o = 3kms — 1 ) at the formation epoch. The 
second and third rows show contour plots of the DF smoothed with observational errors <r b = Skms" 1 , from two realizations with 
the same spiral-structure parameters but different sets of stars. The differences between the second and third rows are a measure of the 
uncertainties in the DFs arising from our numerics. The contours are at 0.2, 0.4, . . . , 1.0 times the maximum value of the DF in each 
panel. 



arises because the typical interval between spiral transients 
in these simulations is 0.5-1 Gyr (Table 0. Thus in most 
of these runs there has been no spiral transient in the last 
~ 0.5 Gyr, and all of the stars formed since the last transient 
still have their birth epicycle energy, which is nearly zero. 
Moreover, these runs feature a relatively small number of 
strong spirals, which hit the disk very infrequently but hard. 
They may therefore overestimate the persistence of moving 
groups. 

To address these concerns, we have run another simu- 
lation sim. In this run, we increase the number of spirals to 
N s = 48 and set e — 0.31. We have also adjusted the central 
times of the spiral transients as follows: (1) We have set the 
central time t s of the most recent transient to be the present 



time to; (2) Because Dehnen's subsample Bl, in which the 
oldest stars have age 0.4 Gyr, already shows two prominent 
moving groups (top left panel of Figure |5J , we have set the 
central time of another transient to to — 0.2 Gyr; (3) The 
other central times are chosen randomly between t s — and 
t 3 — to — 0.2 Gyr. This strategy suppresses the strong mov- 
ing group at the origin and makes the other features in the 
DF more prominent. 

Figure |5] shows the velocity distribution in run sima as 
a function of stellar age. "Moving groups" (i.e. the same 
local concentration of stars in velocity space) can be seen 
among stars of any age, even the oldest stars. Generally with 
increasing age there are more but weaker moving groups, 
however, a strong "moving group" can be present even in a 
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very old sample (see panel 17, 8.5 Gyr). Notice that the same 
"moving group" can be present in stars with a wide range of 
ages; for example, the moving group "B" in panel 13 can be 
traced from panels 11 to 14, a range of at least 1.5 Gyr in 
age. This result is consistent with the observation that the 
same moving grou p is seen in ob servational sub-samples of 
different age (e.g.. |Pehnenlll998l Table 2). 

In Figure ED we show the observed and simulated DFs, 
divided into age bins. The obse rved DF, shown in the top 
row, is taken from lDehnenl jl998T) , The panels represent sub- 
samples determined by color cuts, which should correspond 
to age cuts. The maximum ages in the three subsamples Bl, 
B2, and B3 (the first, second, a nd fourth pan els) are 0.4 
Gyr, 2 Gyr, and 8 Gyr (Table 1 of |Pehnenlll99l) . The sam- 
ple B4GI (fifth panel) is a combination of red (B — V > 0.6) 
main-sequence stars and giant stars, which should have a 
uniform distribution of ages if the star-formation rate has 
been uniform. The third panel, labeled B23, is obtained by 
subtracting the scaled DF of B2 from the DF of B3, in order 
to represent approximately the DF of stars with ages 2-8 
Gyr. 

In the bottom two rows of Figure [5] we show the sim- 
ulated DFs from two runs: run sim, and a modified version 
of sim in which the width ao of the epicycle energy at for- 
mation (eg. I15H is increased from 3kms _1 to Skms" 1 . The 
different panels show the DFs corresponding to the same age 
ranges as in the subsamples of the data. The simulated DFs 
reproduce most of the prominent features of the observed 
DFs: (1) The stars clump together in moving groups. (2) 
The number of moving groups increases as the mean stellar 
age increases, from left to right in the figure. (3) The same 
moving groups are present in stars with a wide range of ages. 
(4) The moving groups appear to be elongated along lines 
of constant v, although this effect is significantly stronger in 
the simulations than the observations. 

Increasing the initial epicycle energy (compare the sec- 
ond and third row of Figure spreads out the moving 
groups somewhat, and enhances their elongation along lines 
of constant v, but they remain prominent; the moments of 
the DF given in Table [5] are almost unchanged. 

In Figure [TU1 we plot the birthplaces of groups "A" and 
"B" shown in panels 2 and 12 of Figure |S| The stars in 
the younger group, A (age 1 Gyr), come from a relatively 
small arc, while the stars in the older moving group B (age 
6 Gyr) were born with a wide range of Galactic azimuths. 
In contrast, Eggen's model for moving groups predicts that 
the stars in a group should be born in a small region and 
thus should have the same age and metallicity. 

We also ran a modified version of 20, which has the 
same parameters as 20 but adjusts the central times of the 
transients as described above for sim. This run gives more 
prominent moving groups than sim, and the groups last for 
a longer time (up to 3.5 Gyr). This result suggests that 
a small number of strong spiral transients tend to produce 
strong moving groups and a lumpy DF, while a large number 
of weak transients tend to produce weak or undetectable 
moving groups and a smooth DF. Therefore, the observed 
moving groups suggest that at most a few tens of spiral 
transients have affected the stars in the solar neighbourhood. 

We also looked for moving groups in run m4, which has 
the same parameters (including the random number seeds) 
as the modified run 20, except that the spirals have four 



arms 2 . We found that substructure in the DF in run m4 
was much weaker than in the modified run 20, presumably 
because the wavelength 2n/k of the arms is smaller so the 
DF is more thoroughly mixed for a given amount of heat- 
ing. The absence of substructure suggests that transient spi- 
rals in the Galaxy may be predominantly two-armed rather 
than four-armed. Further exploration of the effects of the 
spiral-structure parameters on the structure of the local DF 
is obviously worthwhile. 

We have also looked for substructure in the velocity DF 
in runs in which the duration a s of the spiral transients 
was varied (Figure [3J . We found that the strength of the 
substructure declines as a s increases; in other words, spiral 
transients with duration much longer than tend to heat 
the disk without generating moving groups. 

We have also investigated heating by GMCs rather than 
spiral structure, and find that thi s process does no t yield 
significant substructure in the DF llDe Simond EoOO). Thus 
the observed strong substructure suggests that spiral struc- 
ture rather than GMCs is the dominant heating mechanism 
in the solar neighbourhood, a conclusion reached by others 
from independent arguments (see SQ. 



5 INHOMOGENEOUS STAR FORMATION 

So far in this paper we have assumed that the star-formation 
rate is uniform in both space and time, and that both the 
AVR and the structure of the velocity distribution in the so- 
lar neighbourhood are due to scattering by the large-scale 
potential fluctuations associated with transient spiral arms. 
The traditional view, which we explore briefly in this Sec- 
tion, is that the structure in the local velocity distribution 
results from spatial and temporal variations in the star- 
formation rate. Of course, some elements of both pictures 
must be correct, since we know both that spiral arms are 
present and that star formation is concentrated in GMCs. 

The observed AVR requires heating by some 
mechanism — whether spiral arms, GMCs, or one of the 
other mechanisms discussed in ^ However, in this Section 
we ignore the details of this heating, to focus on the 
properties of the irregularities in the velocity distribution 
that are due to inhomogeneous star formation alone. We 
consider the following crude model: 

(i) All stars form in clusters of N stars, of limited spatial 
extent. 



2 It is striking that run m4 yields almost the same velocity disper- 
sion as run 20 (Table|5J. The reason is as follows: consider a simu- 
lation with azimuthal wavenumber m = 2 that yields a power-law 
age- velocity dispersion relation of the form 1291 . <r(t) = (Ct) p . We 
expect that the constant C is proportional to the diffusion coef- 
ficient, and that this in turn is proportional to the mean-square 
fluctuating force F 2 = (fc4> a ) 2 ; thus C = cF 2 where c is a con- 
stant and the age-velocity dispersion relation can be written as 
a(t) = F 2p (ct) p . According to the scaling relations at the end 
of i|3.1l a simulation with azimuthal wavenumber m is equivalent 
to one with m' = fm (here / = 2) if a' = a/f and F' = F/f. 
Thus a'{t) = F 2 P(ct) p /f = fP" 1 {F') 2p (ct) p . Since the exponent 
p ~ 0.5 for run 20, / 2p_1 ~ 1 and we expect the age- velocity dis- 
persion relations for runs 20 and m4 to be nearly the same. 
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Figure 8. The evolution of the velocity distribution as a function of age, from run sima. Each panel is labeled with the age in Gyr. The 
gray scale in each panel is relative to the maximum density in that panel. Notice that the moving group B can be seen at least from 
panels 7 to 13, corresponding to ages of 3.5—6.5 Gyr. 



(ii) The cluster dissociates after the stars form. When the 
cluster dissociates, the DF of its stars is Gaussian in both po- 
sition and velocity, with dispersions a q and a p respectively. 
The stars do not interact with one another once the cluster 
has dissolved; thus, their motion is governed by equations 

Typical values of a q and a p are a few pc or kms" . 

(iii) The heating mechanism does not disrupt or spread 
the phase-space stream of stars from the dissolved cluster 
(any such effects will make many of the conclusions below 
even stronger). We ignore the details of the heating mecha- 
nism, and account for heating by simply assuming that the 
dissociated cluster is born at position (xo,yo), on a non- 
circular orbit with velocities uo and Vo . Thus, the initial DF 
of the dissociated cluster is 



/o(z) 



N 



(2ir<j p a q ) 2 



exp [-i(z - zo)Q(z I 



zo) 



(32) 



where z = (x,y,u,v), T denotes the transpose, and 

q=i I "i :; i 



The covariance matrix is 

<(z T -z r )®(z-z )) = Q- 1 , (34) 
where the average (•) is over the DF /o(z). 
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The trajectories <| 5 } can be rewritten in the form 
jAsiain. Figueras fc Torralll999l) 



s(t) = A(t)z(0), 



(35) 



where t is the time elapsed since the dissociation of the clus- 
ter, and (recall that u — —x) 
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Figure 9. The observed DF as a function of stellar age (Dolmen 1998), compared to simulated DFs from run sima. From left to right, the 
samples represent stars younger than 0.4 Gyr, stars younger than 2 Gyr, stars from 2—8 Gyr in age, stars younger than 8 Gyr, and stars 
of all ages (see text for further details). The bottom row is from a run with the same parameters as sima except that initial dispersion 
erg is 8kms rather than 3kms _1 . 
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Here r = nt. The matrix A(t) has the following properties: 
A(0) = I (by definition), |A(i)| = 1 (by Liouville's theo- 
rem), and K~ x (t) — A(—t) (by time-reversibility). 

Since phase-space density is conserved along the particle 
trajectories, the initial DF is Gaussian (eg. 1321 . and z(t) is a 
linear function of z(0), the DF of the dissociated cluster will 
remain Gaussian for all time. The centroid of the Gaussian 
will be z c (t) = zoA T (t) and the covariance matrix will be 

C(t) = ([z T - x T c {t)] ® [z - z e (i)]> = A(t)Q- 1 A T (t). (37) 

The disrupted cluster is greatly elongated in the az- 
imuthal direction. Therefore to estimate the maximum num- 
ber of stars from a cluster that could be visible in the 
solar neighbourhood, it is useful to compute the linear 
number density at time t, n(y,t) = f f(z,t)dx dv, where 



z = (x, y, u, v) = (x, v). Since the DF is Gaussian, the linear 
density must also be Gaussian, 

N 



n(y,t) 



(27TCT 



2^172" ex P [-(y -ycf/ 2a l] . 



where 



([y - Vc(t)f) = C 22 (t); 



(38) 



(39) 



at la rge times, kJ > 1, we have ijAsiain. Fieueras fc Torral 
1999) 



4A£lt,_ 2 , ,„2^2xl/2 



■ (cr p + 4tt G q 



(40) 



For a typical star- forming cluster, Qa q <C o p (Qu q — 
0.03kms _1 (cr 9 /l pc) compared to typical values a p — 
lkms" 1 ). Thus we may drop the term involving a q , so the 
maximum linear density from a disrupted cluster of age t is 

( ^ 



c(t) 



N 



(2n) 1 / 2 a p t \4:AD. 



(41) 



where the factor in brackets is unity for a fiat rotation curve. 
In the solar neighbourhood, 
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Figure 10. The birthplaces of stars in moving groups "A" and "B" in Figure [S] The star shows the Galactic centre. The stars in a 
moving group are not born at a single location. The locations are obtained from the coordinates of the stars in the sheared sheet by 
equating the polar coordinates R/Rq and <p to 1 + x and y, respectively. 



c (t) = 4 x 10 _5 pc _1 AT 



lkms _1 \ /lOGyr 



(42) 



With the same assumptions, the RMS spatial extent of the 
cluster in the radial direction is 



a x = -Jr [5 sin 2 Kt + (1 — cos fti) 2 ] 



1/2 



(43) 



averaging over the epicycle period we have a x = 
(£) 1/a <V n = 50pc(cr p /lkms- 1 ). 

Now suppose that our survey includes all stars in the 
solar neighbourhood to a distance L. Then if L>a x we can 
detect up to 2n max L stars in the moving group, while if 



L <§C &x we will detect up to (tt/2) 



C L jo x stars in 



the group (these upper limits are achieved if the Sun lies 
exactly on the centroid x c (t)). Extrapolating these formulas 
until they meet, we find that the upper limit to the number 
of detected stars in the moving group is given by 



n mg < 8 x 10 3 iV— -i- 

100 pc <7 p t 

for L > 80 pc(o>/l kms" 1 ), and 



n mg < 1.0 x 10" N 



100 pc 



'lkms^V lOGyr 



(44) 



(45) 



for L < 80pc(cr p /lkms _1 ). The sample of lDehnenl Jl99ST) 
was limited to distances L ~ 100 pc (the median parallax 
error in the Hipparcos catalog is ov = 1 mas and Dehnen's 
sample is restricted to stars with a n /ir < 0.1). If we assume 
that at least n mg = 50 stars are needed to define a moving 
group, then any such group identified in Dehnen's catalog 
must have come from a cluster with initial population 



N > 6 x 10 d 



N > 5 x 10' 



lkms- 1 lOGyr 



if a p < 1.3 km s"\ (46) 



Vlkms- 1 



t 



lOGyr 



if cr p > 1.3 km s" 1 . (47) 



These are conservative lower limits, since (i) we have ne- 
glected the spreading of the stream due to interactions with 
small-scale fluctuations in the gravitational potential, such 
as GMCs; (ii) we have assumed that the radial coordinate of 
the stream of stars is centred precisely on the Sun. 

Infrared observations show that most star formation oc- 
curs in embedded star clusters within GMCs. These clusters 
have masses in the range 5O-1O 3 M0 llLada fc Ladal I2003T) 
and virial dispersions a p ~ 1 kms" 1 and thus are too small 
to contribute detectable moving groups with an age exceed- 
ing ~ 1 Gyr. An alternative is to suppose that the dissoci- 
ated cluster is the GMC itself. In this case a p ~ 5kms~ 
which requires N > 1.3 x 10 5 ; given that the largest GMCs 
have gas masses ~ 10 5 -10 6 Mq we would require extremely 
high star-formation efficiency to produce a sufficient num- 
ber of stars. In other words it is difficult to find a plausible 
astrophysical source for the star-forming regions which are 
supposed to be the progenitors of the old moving groups. 



6 MIGRATION 

ISellwood fc Binnevl <|2002F) discuss the dynamics of radial 
migration due to transient spiral arms, and the relationship 
between radial migration and disk heating. They estimate 
that old stars formed in the solar neighbourhood should now 
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be scattered nearly uniformly throughout the annulus from 
Ro ± 4 kpc. 

We have examined the radial migration or diffusion of 
stars that is induced by spiral-wave scattering, using the 
spiral-arm parameters for our runs 5, 10, 20 and 40. In each 
case we started 100,000 test stars in circular orbits near 
the solar neighbourhood (±0.2 kpc) and followed them for 
10 Gyr. In Figure [Til we plot the RMS change in guiding- 
centre radius and the distribution of stars as functions of 
final epicycle energy. The present in-plane epicycle energy 
of the Sun is only E x = 64(kms~ 1 ) 2 , which is nearly zero 
on the scale of these plots. If these simulations correctly rep- 
resent the properties of the local spiral structure, we expect 
that the radial migration of the Sun since its formation 4.5 
Gyr ago is in the range 0.45 ± 0.25 kpc (pitch angle 10°) to 
1.4 ± 0.9 kpc (pitch angle 20°). The sheared sheet is sym- 
metric in ±x so our simulations cannot address the question 
of whether the migration is likely to be inward or outward; 
ISellwood fc Binnevl J2002D argue that migration is generally 
outward because of the surface-density gradient in the disk. 
These estimates of the migration are consistent w ith the 
estimate made bv lWielen. Fuchs fc DettbarrJ IU996T) . on the 
basis of the Sun's metallicity and the radial metallicity gradi- 
ent in the Galactic disk, that the Sun has migrated outward 
by 1.9 ± 0.9 kpc in the past 4.5 Gyr. 



7 DISCUSSION 

We have explored the hypothesis that transient spiral 
arms are the dominant mechanism that drives the evo- 
lution of the velocity distribut i on in the solar neighbour- 
hood jBarbanis fc Wolti^l967tlSeilwood fc Carlberdll984t 
ICarlberg fc Sellwoodlll985ir The most thorough investiga- 
tions of this mechanism so far have been based on the 
Fokker-Planck equation iJenkins fc Binnevl Il990l : Ijenkinsl 
1992). Our investigation is based on direct numerical in- 
tegration of test-particle orbits in the sheared sheet (eqs. 
|3J. We impose a randomly distributed sequence of trailing 
m — 2, 4 spiral waves with pitch angle a and a strength that 
is Gaussian in time (eq. 12 U , and apply the boundary con- 
ditions that the test particles must be on circular orbits at 
their formation time and in the solar neighbourhood at the 
present time. 

We confirm that transient spiral waves can heat the 
galactic disk in velocity space. The configurations of the spi- 
ral waves in our simulations, such as the number (N„ — 9- 
48), duration (fla s = 0.5-2.5), RMS fluctuation amplitude 
(e rms = 0.16-0.86), and pitch angle (a = 5°-40°) are all 
plausible. Spiral arms can lead to a wide range in the ob- 
served exponent of the age-velocity dispersion relation (eq. 
I29H . p — 0.2-0.76, so this exponent is not a strong discrim- 
inator between different mechanisms of disk heating. The 
variations of the vertex deviation, mean radial velocity, and 
axis ratios of the velocity ellipsoid in the solar neighbour- 
hood are consistent with observations. 

The small-scale structure of the local stellar velocity 
distribution cannot be interpreted in terms of axisymmet- 
ric, steady-state models. Our simulated distribution func- 
tions do not have the Schwarzschild form Q and in addi- 
tion exhibit rich small-scale structure, similar to the "mov- 
ing groups" and "branches" that observers describe in the 



solar neighbourhood distribution function. This result sug- 
gests that moving groups arise primarily from smooth star 
formation in a lumpy potential, in contrast to the tradi- 
tional assumption that they reflect lumpy star formation in 
a smooth potential — although of course elements of both 
pictures must be present to some extent. 

All these findings strongly support the conclusion that 
transient spiral waves are the dominant heating mechanism 
of the disk in the solar neighbourhood. 

Strong substructure in the local velocity distribution is 
most easily produced when the Galaxy experiences a small 
number of short-lived but strong (fractional surface-density 
amplitude e ~ 1) spiral transients, as might be caused by 
minor mergers. A model in which spiral transients are con- 
tinuously present at a lower amplitude could reproduce the 
age- velocity dispersion relation without significant substruc- 
ture. 

We also find some other results: 



(i) The age-velocity dispersion relation depends strongly 
on the stochastic properties of the spiral waves, especially in 
cases where the heating is due mainly to a relatively small 
number of spiral transients (e.g. ~ 10 in runs 20 and 40). 
Therefore, it is difficult, if not impossible, to infer the de- 
tails of the heating mechanism from observations of the age- 
velocity dispersion relation, no matter how accurate. 

(ii) Spiral arms lead to radial migration of stars 
iSellwood fc Binnevl 120021 and references therein). For the 
spiral-wave parameters we used, the RMS radial migration 
of a star like the Sun is a strong function of pitch angle a, 
ranging from 0.45 ± 0.25 kpc for a = 10° to 1.4 ± 0.9 kpc for 
a = 20°. In other words the relation between the amount 
of heating of the solar neighborhood and the amount of ra- 
di al migration is not unique , a conclusion already reached 
bv lSellwood fc Binnevl (2002) using other arguments. 

(iii) Stars in old moving groups did not form at a common 
place and time, and therefore do not necessarily share a 
common age and metallicity. 

(iv) The traditional model, in which structure in the ve- 
locity distribution function is the result of inhomogeneous 
star formation, is difficult to reconcile with at least two ob- 
servations: 1. In this model, all of the stars in a moving group 
should have the same age and metallicity, while in fact the 
stars in prominent moving groups show a wide range of ages. 
2. It is difficult to find a suitable astrophysical source that 
is small enough, cold enough, and produces enough stars to 
be the progenitor of the old moving groups. 

(v) We have also investigated heating by giant molecu- 
lar clouds alone, and find that this process produces much 
weaker substructure in th e distribution function than is ob- 
served (iDe Simond 12000): thus giant molecular clouds are 
unlikely to be the dominant heating mechanism in the solar 
neighbourhood. 



Although heating by spiral waves explains most of the 
properties of the solar neighbourhood velocity distribution, 
there is still room for other heating mechanisms, such as the 
Galactic bar, giant molecular clouds, and halo substructure 
or minor mergers. 
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Figure 11. (top) The mean of |Ai?|, the absolute value of the change in guiding-centre radius, as a function of final epicycle energy 
Ex = ^K 2 a 2 , where a is the epicycle amplitude (eq. I1UI . The stars are initially in circular orbits with radius Rq. The four panels 
correspond to scattering by spiral arms with the parameters of runs 5, 10, 20 and 40 (Table 0. The error bars represent the standard 
deviation in the distribution of |AH| at a given radial energy, (bottom) The distribution of stars as a function of epicycle energy. Dashed 
lines are the Schwarzschild distribution with the same velocity dispersion. 
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